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Abstract. We discuss recent advances made in modelling the complex 
magnetohydrodynamics of the Sun using our anelastic spherical harmon- 
ics (ASH) code. We have conducted extensive 3-D simulations of com- 
pressible convection in rotating spherical shells with and without mag- 
netic fields, to study the coupling between global-scale convection and 
rotation in seeking to understand how the solar differential rotation is 
established and maintained. Such simulations capable of studying fairly 
turbulent convection have been enabled by massively parallel supercom- 
puters. The resulting convection within domains that capture a good 
fraction of the bulk of the convection zone is highly time dependent and 
intricate, and is dominated by intermittent upflows and networks of strong 
downflows. A high degree of coherent structures involving downflowing 
plumes can be embedded in otherwise chaotic flow fields. These vortical 
structures play a significant role in yielding Reynolds stresses that serve 
to redistribute angular momentum, leading to differential rotation pro- 
files with pole to equator contrasts of about 30% in angular velocity Q 
and some constancy along radial line at mid latitudes, thereby making 
good contact with deductions from helioseismology. When a magnetic 
field is introduced, a dynamo regime can be found that does not destroy 
the strong differential rotation achieved in pure hydrodynamics cases. 
The magnetic fields are found to concentrate around the downflowing 
networks and to have significant north-south asymmetry and helicity. 



1. Observational Challenges 



The surface layers of the Sun have long been known to exhibit complex con- 
vection and magnetism involving a very broad range of spatial and temporal 
scales (Stix 2002). Helioseismology is now permitting novel views of the interior 
structure and dynamics within our nearest star (Gough & Toomre 1991). Using 
millions of acoustic modes, it is possible to probe as a function of radius and lat- 
itude the solar sound speed c, density p and angular velocity £1. Figure 1 shows 
the solar internal rotation profile as a contour plot and as six latitudinal cuts 
(Schou et al. 1998, Howe et al. 2000). Although there are prominent variations 
of O with latitude, with the equator rotating considerably faster than the high 
latitudes, it is quite striking that Q is largely constant at mid latitude along 
radial lines and imprints itself all the way down to the base of the convective 
zone. There a tachocline of strong velocity shear joins the nearly solid body 
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rotation of the radiative interior with the differentially rotating convective zone. 
A distinctive near-surface shear layer is also evident. Such angular velocity pat- 
terns were not anticipated and are at variance with earlier models of differential 
rotation within convecting rotating shells (Glatzmaier & Gilman 1982, Gilman 
& Miller 1986, Glatzmaier 1987). In those models O was mostly constant along 
cylinders aligned with the rotation axis (e.g. Taylor columns). 




Figure 1. (a) Angular velocity profile Q/2n in latitude and radius 
as deduced through inversion of acoustic mode frequency splittings 
from the SOI-MDI helioseismic instrument; the equator here coincides 
with the horizontal axis [adapted from Schou et al. 1998]. (b) Time- 
averaged rotation rates from five years of GONG helioseismic data, 
plotted against proportional radius at different latitudes, with rapid 
rotation at the equator and slower rotation at high latitudes. The zone 
covered by our computational domain is indicated (grey area) [adapted 
from Howe et al. 2000]. 

Intimately related to the dynamics of the solar turbulent convection zone 
is its magnetic activity, variously involving sunspots, prominences and CME's, 
along with its overall 22-year cycle. How such a turbulent and complex system as 
the Sun can exhibit order amidst what is seemingly chaos is a most challenging 
question. It is generally thought that the solar magnetic dynamo operates at 
two differing ranges of spatial and temporal scales (Cattaneo & Hughes 2001). 
The global dynamo yielding the regular 22-year cycle and butterfly diagrams 
for sunspot emergence is likely to be seated within the tachocline at the base 
of the convection zone (Parker 1993). The origin of the rapidly varying and 
smaller scale magnetism is probably due to local dynamo action achieved by the 
intensely turbulent convection. We would like here to cast new light on some 
aspects of this complex magnetohydrodynamical (MHD) system. We believe 
that 3-D MHD numerical simulations are essential to pursue questions of solar 
magnetism. Given the large range of temporal and spatial scales involved in 
the solar convective envelope, one has to choose between either a local high 
resolution domain or a global and somewhat less resolved spherical domain. The 
advantage of the former is its ability to resolve highly turbulent flows (Brummell 
et al. 1998, 2002; Stein & Nordlund 1998). The global approach, by fixing the 
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largest scales, is more limited in the turbulence levels that can be resolved, 
yet takes into account the correct geometry and its topological implications for 
mean flows. We will focus our attention here on the establishment of the solar 
differential rotation and meridional circulations, both of which are the purview 
of global models. We shall present recent results from our simulations of solar 
convection within full 3-D spherical shells, discussing also the angular velocity 
Q profiles that can be achieved in the bulk of the convection zone and the level 
of dynamo induced magnetism that can be sustained there. 



2. Our Numerical Approach 

The ASH code solves the 3-D MHD anelastic equations of motion in a rotating 
spherical shell geometry (Clune et al. 1999, Miesch et al. 2000). These equations 
are fully nonlinear in velocity and magnetic field variables; the thermodynamic 
variables are separated with respect to a spherically symmetric and evolving 
mean state having a density p, pressure P, temperature T and specific entropy 
S, and fluctuations about this mean state, namely p, P, T, S: 



V-(pv) = 0, (1) 

/<9v \ 1 

pi — + (v- V)v + 2r2 Q x vj = —VP + pg + — (V x B) x B 

- V-X>-[VP-pg], (2) 
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pf— = V-[K r pc p V(f + T) + K pfV(S + S)} 

-pTvV(S + S) + ^(VxB) 2 + 2pu [etjeij - 1/3(V • v) 2 ] , (3) 
9B 

— = V x (v x B) - Vx (r?V x B), (4) 

where v = (v r ,V0,v^) is the local velocity in spherical coordinates in the frame 
rotating at constant angular velocity fi , g is the gravitational acceleration, 
B is the magnetic field, c p is the specific heat at constant pressure, n r is the 
radiative diffusivity, rj is the effective magnetic diffusivity, and T> is the viscous 
stress tensor, involving the components 

V i j = -2p v [e i j-l/S(V-v)8 i j], (5) 

where is the strain rate tensor, and v and re are effective eddy diffusivities. 
To complete the set of equations, we use the linearized equation of state 

p P f 1 P cp' U 
where 7 is the adiabatic exponent, and assume the ideal gas law 

P = Tlpf (7) 

where 1Z is the gas constant. 
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The velocity and magnetic fields and the thermodynamic variables are ex- 
panded in spherical harmonics for their horizontal structure and in Chebyshev 
polynomials for their radial structure. This approach has the advantage that 
the spatial resolution is uniform everywhere on a sphere when a complete set of 
spherical harmonics is used up to some maximum in degree I (retaining all az- 
imuthal orders m). The anelastic approximation captures the effects of density 
stratification without having to resolve sound waves which would severely limit 
the time steps. We use a toroidal and poloidal decomposition that enforces the 
mass flux and the magnetic field to remain divergence free. 

The model is a highly simplified description of the solar convection zone: 
solar values are taken for the heat flux, rotation rate, mass and radius, and a 
perfect gas is assumed since the upper boundary of the shell lies below the H 
and He ionization zone. The computational domain extends from 0.72 R® to 
0.96 or 0.98 R®, thereby concentrating on the bulk of the unstable zone and 
here not dealing with penetration into the radiative interior. The effects of the 
steep entropy gradient close to the surface has been softened by introducing a 
subgrid scale (SGS) transport of heat to account for the unresolved motions, 
and enhanced eddy diffusivities are used in these large eddy simulations (LES). 
The typical density difference across the shell in radius is about 30. 

In order to get reliable statistics, these numerical experiments need to be 
integrated over long periods in physical time. Within our 3-D convection sim- 
ulations the dynamical time is of the order of 30 days (the time roughly taken 
by a fluid element to travel across the convection zone). This implies that time 
averages of relevant quantities such as fi have to be performed over no less than 
300 days when feasible. We here present one of our most turbulent cases, namely 
case E computed with a Prandtl number P r = v j k of 0.25 and a rms Reynolds 
number R e = v rms D/v of ~ 650, where D is the thickness of the shell. This 
case has been restarted from a promising laminar case AB (Brun & Toomre 
2002) which had been evolved over 200,000 time steps for a total physical time 
of about 6200 days starting from quiescent initial conditions. To achieve more 
turbulent states such as in case E, we have lowered the viscous and thermal 
effective diffusivites v and k in case AB gradually in a sequence of steps, wait- 
ing for the simulation to relax between each step. For case E we present the 
latest 85 days of the simulation over a total of 9,200 physical days of evolution 
(corresponding to times steps 430,000 to 480,000). At this level of turbulence, 
numerical accuracy requires a resolution of £ max = 680 and a physical time step 
of 150 s. These numbers illustrate the difficulty of conducting such 3-D tur- 
bulent simulations. Although our ASH code is sustaining 300 Mflops/cpus or 
20% of peak performance, thus performing very well on massively parallel su- 
percomputers such as the IBM-SP3 at SDSC, such a single simulation requires 
of order 80,000 cpu hours. Typical runs use 256 processors, which is equivalent 
to 76 Gflop/s. We now describe in details some results from our 3-D convection 
simulations. 



3. Turbulent Convection under the Influence of Rotation 

Figure 2 displays the evolution of radial velocity in case E over 10 days in time 
near the top of the domain. Dark tones represent downflows and lighter ones 
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upflows. The vantage point is in the uniformly rotating frame used in our simu- 
lations. The convection patterns are intricate and highly time dependent. Some 
of the pattern evolution is related to the advection by prograde zonal flows near 
the equator and retrograde ones at higher latitudes associated with the differ- 
ential rotation driven by the convection relative to this frame. Shearing and 
cleaving of the convective cells is evident, as well as distortions in the downflow 
lanes. There is an asymmetry between the broad upflows centered in each con- 
verting cell and the narrow fast downflows at their periphery. This leads to a 
downward transport of kinetic energy. The strong correlations between warm 
upward motions and cool downward motions are essential in transporting the 
heat outward. 




Figure 2. Evolution in the convection over 10 days, showing the ra- 
dial velocity in case E near the top (0.97 Rq) of the spherical domain. 
The time interval between each successive image is about 5 days. Down- 
flows appear dark. The dotted circle is located at the solar radius Rq 
and the equator is indicated by the dashed curve. 



Pronounced vortical structures are evident at the interstices of the down- 
flows network. They are counterclockwise in the northern hemisphere and clock- 
wise in the southern one, i.e cyclonic. The strongest of these vortex tubes or 
'plumes' extend through the whole domain depth. These plumes represent co- 
herent structures that are surrounded by more chaotic flows. They tend to align 
with the rotation axis and to be tilted away from the meridional planes, leading 
to Reynolds stresses that are crucial ingredients in redistributing the angular 
momentum within the shell. 

3.1. Making Good Contact with the Solar Differential Rotation 

The differential rotation profile in latitude and radius associated with the vig- 
orous convection of case E is shown in Figure 3. For simplicity, we have con- 
verted the mean longitudinal velocity va, into a sidereal angular velocity $7, us- 
ing £l /2ir = 414 nHz (or 28 days) as the reference frame rotation rate. In 
the contour plot, the near polar regions have been omitted due to the difficulty 
of forming stable averages there, since the averaging domain is small but the 
temporal variations large. Case E exhibits a fast (prograde) equatorial region 
and slow (retrograde) high latitude region. This is due to correlations in the 
velocity components leading to significant Reynolds stresses. These Reynolds 
stresses are intimately linked to the influence of Coriolis forces acting upon the 
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Figure 3. Temporal and longitudinal averages in cases E of the an- 
gular velocity profiles formed over an interval of 85 days. This case 
exhibits a prograde equatorial rotation and a strong contrast AQ from 
equator to pole, as well as possessing a high latitude region of par- 
ticularly slow rotation. In the right panel a sense of the asymmetry 
present in the solution can be assessed in these radial cuts at indicated 
latitudes. 

convecting motions and to the presence of plumes tilted both away from the 
local radial direction and out of the meridional plane. Such correlations have 
been identified in local high resolution Cartesian domains as well (Brummell 
et al. 1998). These lead to an equatorward transport of angular momentum, 
resulting in the slowing down of the high latitude regions and speeding up of 
the equatorial zone. At low latitudes there is some alignment of 12 along the 
rotation axis. At mid latitudes, the angular velocity is nearly constant along ra- 
dial lines, in good agreement with helioseismic deductions (cf. Fig. 1). Further, 
case E exhibits a monotonic decrease of with latitude, a property that has 
been difficult to achieve in 3-D spherical convection calculations. Indeed, most 
other cases have their equator to pole contrast Af2 confined to mid latitudes (i.e 
where the inner tangent cylinder cuts through the outer shell at ~ 42°). The 
differential rotation contrast between the equator and 60° in case E is 110 nHz 
(or 26% relative to the frame of reference), thus being very close to the 92 nHz 
(or 22%) variation observed in the Sun. Since our progenitor case AB shared 
this attribute as well, it is comforting that we can retain this solar-like property 
in a significantly more turbulent and complex case such as case E. Many of the 
more complex cases that we have computed previously had a tendency to lose 
some of their latitudinal contrast in £1, which also became more nearly constant 
along cylinders aligned with the rotation axis (Brun & Toomre 2002). A sense of 
the asymmetry present in case E can be assessed both in the contour plot and in 
the latitudinal cuts (right panel of Fig. 3) , where we have plotted Q in the north 
(dotted) and south (dashed) hemispheres along with their mean. The convection 
itself exhibits some asymmetry between the two hemispheres (cf. Fig. 2), and 
so it is not surprising that the mean flows driven by the convection do the same. 



APS Conf. Ser. Style 



7 



These asymmetries are expected to diminish over a longer temporal average. 
Mean field models of the solar differential rotation (Kichatinov & Riidiger 1995, 
Durney 1999) have advocated that a thermal wind balance (involving pole to 
equator temperature contrasts) could be the cause of the non-cylindrical profile 
in fi. This could come about through the baroclinic nature of the convecting 
motions yielding some latitudinal heat fluxes, resulting in the breakdown of the 
Taylor-Proudman theorem (Pedlosky 1987). Although it is indeed true that case 
E exhibits latitudinal variation of entropy and temperature fluctuations relative 
to the mean, these are not the most dominant players everywhere in the shell. 
A temperature contrast of few degree K seems compatible with a AQ/£l of 
~ 30%. However, we find that the Reynolds stresses are the main agents re- 
sponsible for the equatorial acceleration achieved in our simulations, and thus 
the solar differential rotation is dynamical in origin. 

3.2. Pending Issues Concerning the Meridional Circulation 

The meridional circulation associated with the vigorous convection in case E 
is maintained variously by Coriolis forces acting on the differential rotation, 
by buoyancy forces, by Reynolds stresses and by pressure gradients, and thus 
can be thought as a small departure from geostrophic balance. The merid- 
ional circulation exhibits a multi-cell structure both in latitude and radius, and 
given the competing processes for its origin, it is not straightforward to predict. 
Typical amplitudes for the velocity are of order 25 m/s, comparable to local 
helioseismic deductions (Haber et al. 2002). The flow is directed poleward at 
low latitudes, with return flow deeper down. The temporal fluctuations in the 
meridional circulation are large and thus stable time averages are only attained 
by sampling many rotations. The kinetic energy in the differential rotation and 
in the convective motions are two orders of magnitude higher than that in the 
meridional circulation (Brun & Toomre 2002). As a result, small fluctuations 
in the convective motions and differential rotation can lead to major variations 
in the circulation. Some of the helioseismic inferences suggest the presence of 
single cell circulations, which are at odds with our multi-cell patterns. However 
these inferences vary from year to year, and there is recent evidence for double- 
cell structure in the circulations observable in the near-surface shear layer, but 
only in the northern hemisphere as the current solar cycle advances (Haber et 
al. 2002). From a careful analysis of the angular momentum transport in our 
shell we have deduced that the slow pole behavior seen in case E seems to come 
about from a relatively weak meridional circulation at high latitudes. Case E 
shares with case AB the property of relatively mild meridional circulation at the 
higher latitudes. This permits a more efficient extraction of angular momentum 
by the Reynolds stresses from the high latitudes toward the equator in yielding 
the interesting differential rotation profile that is achieved. 

4. Adding Magnetic Fields to the Solar Cauldron 

We now turn to consider the influence of magnetic fields both upon the con- 
vection in our deep shell and upon the angular velocity profiles that can be 
maintained. Early attemps to explain the 22-year solar cycle considered the 
possibility that the solar dynamo operated within the bulk of the convective 
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envelope (Gilman 1983, Glaztmaier 1985), but such approaches failed because 
strong magnetic fields could not be stored efficiently within the unstable stratifi- 
cation of the convection zone. The dynamo periods were too short (of the order 
of 1 year) and the poloidal fields were found to propagate poleward, at variance 
with observations. More recently, Parker (1993) has proposed an interfacial 
dynamo model seated in the stable tachocline (see also Riidiger & Branden- 
burg 1995). Magnetic field is still generated within the bulk of the convection 
zone, but is pumped downward into the stable layer via overshooting turbulent 
plumes, to be stretched there into large-scale toroidal structures by the strong 
shear of the tachocline. When the amplification of the toroidal magnetic field is 
great enough, the structures (called magnetic flux tubes in mean field models) 
become magnetically buoyant and rise upward through the convective envelope. 
The strongest of those structures emerge in the photosphere as bipolar mag- 
netic arcades, whereas the weaker ones are recycled within the convective zone. 
This leads to the crucial natural cycle of poloidal to toroidal interchange, i.e. 
Bp i — > B tor — * B po i. However, many aspects of these essential 'dynamo build- 
ing blocks' remain to be demonstrated through nonlinear 3-D calculations. At 
present it is not feasible to simulate self consistently all the processes operating 
together, and thus one needs to concentrate on individual components (magnetic 
generation, pumping, shearing, buoyant rising). One important ingredient in the 
interfacial dynamo scenario is the ability of the convective motions to generate 
and sustain magnetic fields in the bulk of the zone. We thus have evaluated some 
conditions in 3-D compressible convection for which such a dynamo threshold 
can be realized. We further wished to identify the maximum nonlinear ampli- 
fication of the magnetic field that can be sustained by the convective motions 
without destroying the strong angular velocity contrasts previously attained. 
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Figure 4. Kinetic energy (KE) and magnetic energy (ME) for cases 
Ml, M2 and M3, involving in turn a magnetic rms Reynolds number 
R em of 250 (dashed dot line), 300 (solid) and 500 (dashed). 
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4.1. Nonlinear Dynamo Threshold and MHD Turbulence 

We have conducted three MHD simulations (named respectively Ml, M2 and 
M3) started from a solution slightly more turbulent than case AB but exhibiting 
a similar angular velocity profile £1. We then introduced a small seed dipolar 
magnetic field and let the simulations proceed. Because we had to solve two extra 
equations (i.e the poloidal and toroidal components of the induction equation 
[4]) and integrate these solutions over several ohmic diffusion time, it was not 
feasible to use our more turbulent case E that required higher spatial resolution. 
Figure 4 shows the magnetic energy evolution for three values of the magnetic 
diffusivity rj (i.e. 2 x 10 12 , 1.6 x 10 12 , 10 12 cm 2 /s). We note that over more 
than 3000 days (corresponding to several ohmic decay times) the two lowest 
diffusive cases M2 and M3 exhibit a sustained magnetic energy (ME), the levels 
of which depend on r\. The other case Ml is clearly decaying, since the rate of 
generation of magnetic fields could not compensate for the rate of destruction by 
ohmic diffusion. The dynamo threshold seems to be around a magnetic Reynolds 
number R em = v vms D/rj of ~ 300 or r] ~ 1.6 x 10 12 . This is about 25% higher 
than in a progenitor incompressible simulation (Gilman 1983). 




Figure 5. Snapshots of the radial component of the velocity (left) 
and magnetic (right) fields in case M3 in the middle of the domain 
(0.86 i?o). Dark tones represent downflow (or negative polarity). The 
velocity and magnetic field peak amplitudes are about 100 m/s and a 
few thousand Gauss respectively. The outer dotted circle corresponds 
to solar radius Rq and the dashed curve indicates the equator. 

Further, we note that the kinetic energy (KE) in model M3 has been re- 
duced by about 40% compared to its initial value. In this case ME has grown 
to reach a value of 7% of KE. The increase of the magnetic energy started to 
influence the total amount of kinetic energy contained in the shell when ME 
reached roughly 0.5% of KE after about 600 days of evolution. The early ex- 
ponential growth of ME in case M3 extends to about 600 days from the start, 
after which the nonlinear feedback of the Lorentz forces on the flow begins to be 
felt. For case M2, ME is still small enough (i.e < 0.1%) even after 4000 days for 
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the convective motions to only be mildly affected by the Lorentz forces. This is 
most clearly seen in comparing the kinetic energy time trace for cases M2 and 
M3. Figure 5 displays snapshots of the radial component of both the velocity 
and the magnetic fields for case M3 in the middle of the domain. The magnetic 
field is mainly concentrated in the downflow lanes, having been swept away from 
the center of the convective cells by the broad upflows. Both polarities coexist 
at the downflow network interstices. Clearly the magnetic field has a finer and 
more intricate structures, exhibiting many swirls. This is due to our choice of 
magnetic diffusivity being the smaller, noting that the magnetic Prandtl num- 
ber P rm = v jr] is 4 in this solution. The magnetic energy of the toroidal field 
within the bulk of the convective zone is roughly an order of magnitude stronger 
than that of the poloidal field. Figure 6 displays a 3-D rendering of the toroidal 
field. Substantial magnetic helicity is present, involving complex winding of the 
toroidal magnetic fields along their length, with both polarities interchanging 
their position into structures not unlike cables. The toroidal magnetic fields 
also possess the greater spatial scales, having been stretched by the gradients in 
angular velocity. Some features could resemble magnetic flux tubes, although 
they are short lived. There is a clear north-south asymmetry in both the toroidal 
and poloidal magnetic field topology. 




Figure 6. Snapshot of the toroidal component of the magnetic field 
for case M3 appearing here as 3-D torus centered around the equator 
and slightly tilted forward. 

4.2. Maintaining a Strong Solar-like Differential Rotation 

With fairly strong magnetic fields sustained within the bulk of the convection 
zone in case M3, it is to be expected that the differential rotation f2 will respond 
to the feedback from the Lorentz forces. Figure 7 shows the time averaged 
angular velocity achieved in case M3. As for the convective motions, the main 
effect of the Lorentz forces is to extract energy from the kinetic energy stored in 
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the differential rotation. The reduction of KE contained in the angular velocity 
is of the same order as the decrease seen in the total KE, i.e. 40%. As a 
consequence the angular velocity contrast AQ from 60° to the equator drops by 
~ 30% in case M3, going from 140 nHz (or 34% compared to the reference frame 
fi ) in the hydrodynamic case to 100 nHz (or 24%). Nevertheless, the angular 
velocity in case M3 remains in good agreement with the solar profile, both in 
amplitude and profile, as can been seen by comparing Fig. 1 and 7. The source 
of the reduction of the latitudinal contrast of f2 can be attributed to the poleward 
transport of angular momentum by the Maxwell stresses (the mean magnetic 
fields having a negligable contribution). Now the Reynolds stresses again need 
to balance the angular momentum transport by the meridional circulation, the 
viscous diffusion and the Maxwell stresses. This leads to a less efficient speeding 
up of the equatorial regions. Since ME is only 7% of KE in case M3, the Maxwell 
stresses are not yet the main players in redistributing the angular momentum. 
We have found that a value of ME above about 20% of KE leads to a significant 
magnetic braking effect on the differential rotation. Had the simulation been 
restarted with a stronger initial magnetic field Bo, AQ could drop by 90% in 
less than a few hundreds days, thus being at variance with helioseismic findings. 
By letting the convective motions gradually adapt themself to the nonlinear 
feedback of the Lorentz forces one seems to limit the equipartitioning of ME 
and KE and thus succesfully retain the strong differential rotation seen in pure 
hydrodynamical cases. 




Figure 7. Temporal and longitudinal averages of the angular veloc- 
ity profiles achieved in case M3 over an interval of 80 days. This case 
exhibits a prograde equatorial rotation and a strong contrast A£l from 
equator to pole, as well as possessing a high latitude region of par- 
ticularly slow rotation. In the right panel a sense of the asymmetry 
present in the solution can be assessed in these radial cuts at indicated 
latitudes. 

Our 3-D simulations of convection in deep spherical shells, achieved through 
use of massively parallel supercomputers, are helping to show how the strong 
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differential rotation present in the Sun may be maintained through fairly com- 
plex redistribution of angular momentum by the turbulent compressible flows. 
We have also begun to study the interaction of convection and rotation with 
seed magnetic fields in such shells, thereby identifying some parameter ranges 
in which sustained magnetic dynamo action can be realized without unduly 
reducing the angular velocity contrasts maintained by the convection. 
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